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ABSTRACT 

Image  registration  is  an  essential  step  in  many  image  processing 
applications  that  need  visual  information  from  multiple  images  for 
comparison,  integration  or  analysis.  Recently  researchers  have  intro¬ 
duced  image  registration  techniques  using  the  log-polar  transform 
(LPT)  for  its  rotation  and  scale  invariant  properties.  However,  the 
accuracy  of  the  approach  limits  by  the  number  of  samples  used  in 
the  mapping  process,  which  affects  directly  the  computational  cost. 
Motivated  by  the  success  of  LPT  based  approach  and  its  limita¬ 
tion,  we  propose  a  novel  Projective  Polar  Transform  (PPT)  based 
image  registration  method  that  is  robust  to  translation,  scale,  and 
rotation  and  yields  high  accuracy  while  requires  low  computational 
cost.  Unlike  LPT  that  2D  interpolation  is  needed  in  the  mapping 
process,  our  method  uses  one-to-one  mapping  mechanism  that  di¬ 
rectly  arranges  image  from  Cartesian  to  Polar  coordinate  according 
to  pre-computed  PPT  map.  An  innovative  projection  mechanism  is 
proposed  to  reduce  the  image  from  two  to  one  dimensional  vectors. 
With  the  proposed  approach,  the  mapping  procedure  is  accelerated 
and  the  matching  process  can  be  performed  in  ID.  This  results  in 
great  reduction  of  the  computational  cost. 

Index  Terms —  Aerial  image  registration,  projective  polar  trans¬ 
form. 

1.  INTRODUCTION 

Image  registration  is  a  process  of  aligning  two  images  that  share 
common  visual  information  such  as  images  of  the  same  object  or 
images  of  the  same  scene  taken  at  different  geometric  viewpoints, 
different  time,  or  by  different  image  sensors.  Image  registration 
is  an  essential  step  in  many  image  processing  applications  that  in¬ 
volve  multiple  images  for  comparison,  integration  or  analysis  such 
as  image  fusion,  image  mosaics,  image  or  scene  change  detection, 
and  medical  imaging.  The  main  objective  of  image  registration  is 
to  find  the  geometric  transformations  of  the  model  image ,  /M,  in 
the  target  image ,  IT ,  where  IT(x,y)  =  T{IM (x' ,y')}  and  T 
is  a  two-dimensional  geometric  transformation  that  associates  the 
(x',  y')  coordinates  in  IM  with  the  (x,  y)  coordinates  in  IT .  These 
two-dimensional  geometric  transformations  include  scale,  rotation, 
and  translation  in  the  Cartesian  coordinates.  Many  works  have  been 
done  in  this  area  in  the  past  20  years  [1],  which  can  be  catego¬ 
rized  into  two  major  groups:  the  feature-based  approach  and  the 
area-based  approach.  The  feature-based  approach  uses  only  the 
correspondence  between  the  features  in  the  two  images  for  regis¬ 
tration.  The  features  can  be  color  gradient,  edges,  geometric  shape 
and  contour,  image  skeleton,  or  feature  points.  Since  only  the  fea¬ 
tures  are  involved  in  the  registration,  the  feature -based  approach  has 


advantages  in  registering  images  that  are  subjected  to  alteration  or 
occlusion.  However,  the  use  of  the  feature-based  approach  is  rec¬ 
ommended  only  when  the  images  contain  enough  distinctive  fea¬ 
tures  [2].  As  a  result,  for  some  applications  such  as  aerial  imaging, 
in  which  features  are  difficult  to  be  distinguished  from  one  another 
(for  example,  shapes  of  most  buildings  from  the  bird-eye  view  are 
almost  identical),  the  feature-based  approach  may  not  perform  effec¬ 
tively.  This  problem  can  be  overcome  by  the  area-based  approach,  in 
which  part  of  the  images  are  used  for  comparison  using  correlation¬ 
like  approaches  such  as  cross-correlation  and  Fourier-based  phase 
correlation. 

Innovative  area-based  approach  using  Log-Polar  transform 
(LPT)  method  have  been  proposed  [3-7]  and  proven  to  be  effective 
and  robust  to  changes  in  translation,  scale,  and  rotation.  However, 
there  is  limitation  to  LPT  based  approach.  The  accuracy  of  the 
registration  depends  directly  on  the  number  of  samples  used  in  the 
mapping  process.  This  problem  can  be  easily  illustrated,  for  ex¬ 
ample,  if  36  samples  in  the  angular  direction  are  used,  the  level  of 
accuracy  of  the  registration  can  not  exceed  10  degrees  resolution. 
For  aerial  image  registration,  in  which  scale  and  rotation  parameters 
between  two  images  may  be  very  small,  large  number  of  samples  in 
both  radius  and  angular  directions  are  required.  This  increases  the 
computational  load  tremendously  from  both  the  mapping  process, 
which  require  2D  interpolation  of  image  pixels  in  the  Cartesian,  and 
the  matching  process  that  involves  computation  of  2D  correlation. 

Motivated  by  the  success  of  LPT  based  approach  and  its  limi¬ 
tation,  we  propose  a  novel  Projective  Polar  Transform  (PPT)  based 
image  registration  method  that  is  robust  to  translation,  scale,  and 
rotation  and  yields  high  accuracy  while  requires  low  computational 
cost.  Unlike  LPT  that  2D  interpolation  is  needed  in  the  mapping 
process,  our  method  uses  one-to-one  mapping  mechanism  that  di¬ 
rectly  arranges  image  from  Cartesian  to  Polar  coordinate  according 
to  pre-computed  PPT  map.  An  innovative  projection  mechanism  is 
proposed  to  reduce  the  image  from  two  to  one  dimensional  vectors. 
With  the  proposed  approach,  the  mapping  procedure  is  accelerated 
and  the  matching  process  can  be  performed  in  ID.  This  results  in 
great  reduction  of  the  computational  cost. 

The  paper  is  organized  as  follows.  In  section  2,  we  discuss  our 
proposed  registration  approach  in  detail.  Experimental  results  are 
shown  in  section  3.  Finally,  our  work  is  summarized  in  section  4. 

2.  PROPOSED  ALGORITHM 
2.1.  Projective  Polar  Transforms 

Inspired  by  the  scale  and  rotation  invariance  properties  of  LPT  based 
approach,  we  propose  a  novel  image  transformation  scheme,  PPT, 
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Fig.  1.  PPT  map  (a)  Distance  parameters  for  image  of  size  10  x  10 
pixels  and  (xc,  yc)  =  (7,  7),  (b)  Graphic  demonstration  of  PPT  map 
for  image  with  size  513  x  513  pixels  and  (xc,  yc)  =  (256,  256) 

that  is  robust  to  translation,  scale,  and  rotation  and  yields  high  accu¬ 
racy  while  requires  low  computational  cost. 

2.1.1.  PPT  mapping 

To  transform  a  circular  area  with  the  size  Rmax  inside  image  I(x,y) 
in  the  Cartesian  coordinate  to  projective  polar,  we  first  need  to  create 
a  PPT  map.  We  define  a  distance  parameter  -D(®c,i/c)  ( x ,  y)  as 

D(xc,yc)(x,y)  =  max{n  G  Z+|n  <  y/(x-  *c)2  +  {y  -  t/c)2} 

(1) 

where  (xc,  yc)  is  the  image  pixel  in  Cartesian  selected  to  be  the  cen¬ 
ter  of  the  transformation.  Fig.  1(a)  shows  example  of  distance  pa¬ 
rameters  of  image  with  the  size  10x10  pixels  where  the  center  point 
(xc,  yc)  =  (7,  7).  In  each  pixel  of  the  image,  the  calculated  distance 
parameter  is  filled.  Fig.  1(b)  is  the  graphic  demonstration  of  PPT 
map  for  image  with  the  size  513  x  513  pixels  where  the  center  of  the 
transformation  (xc,  yc)  =  (256,  256).  As  shown  in  Fig.  1(b),  PPT 
map  is  an  approximation  of  Polar  transform  mapping. 

To  transform  image,  we  directly  map  the  image  pixels  in 
the  Cartesian  coordinate  I(x,y)  that  have  distance  parameter 
D(Xc,yc)(x,y)  —  i  to  the  transformed  coordinate  IP(i,9).  To 
approximate  sampling  in  both  radius  and  angular  directions,  for 
each  iteration  i,  the  mapping  begins  with  I(xc,  yc  +  i)  and  continue 
the  search  for  image  pixels  with  D(Xc^cfx,  y)  =  i  pixel-by-pixel 
in  the  counter-clockwise  direction.  The  process  is  repeated  for 
i  =  1, ..,  Rmax •  Example  of  the  transformed  Lena  image  is  shown 
in  Fig.  2.  Since  PPT  map  is  static  and  does  not  required  to  be  com¬ 
puted  for  every  transformation,  and  the  mapping  process  is  a  simple 
one-to-one  allocation  from  Cartesian  to  PPT  images  and  does  not 
require  interpolation,  the  computational  cost  in  this  step  is  very  low 
compared  to  that  of  LPT. 

2.1.2.  Projection  transform 

As  shown  in  Figs.  2(b),  the  result  from  the  transformation  is  a  series 
of  sample  bins  arranged  in  the  step-like  manner  which  do  not  show 
coordinate  shifts  for  scaled  and  rotated  image  as  in  LPT  any  longer. 
To  maintain  the  advantages  of  scale  and  rotation  invariance  in  LPT, 
we  use  an  innovative  projection  transform  method  which  projects 
the  two-dimensional  image  on  the  radius  and  angular  coordinates, 
respectively.  From  the  two  projections,  we  can  accurately  calculate 
the  scale  and  rotation  parameters,  for  which  the  details  will  be  elab¬ 
orated  in  Section  2.2.  For  now  we  define  the  projection  transform 
for  the  APT  transformed  image. 


Given  a  transformed  image  IP(r ,  6)  that  consists  of  nr  sample 
bins  (nr  =  Rmax),  in  which  each  bin  has  the  length  of  nei  for 
i  =  1, ...,  nr.  We  denote  rio ,  3ft  and  0  as  the  number  of  samples  in 
the  angular  direction  at  n  =  Rmax,  the  projection  on  the  radius  co¬ 
ordinate,  and  the  projection  on  the  angular  coordinate,  respectively. 
The  mathematical  expressions  of  3ft  and  ©  are  as  follow: 

nei 

=  (2) 
3  = 1 

©O')  =  (3) 

i—  1 

i  e  [1,  ...,nr],  j  €  [1, rio],  = 

Vij  =  n*0  “  !)  “ floor[Sli(j  -1)],  riij  =  l-rilj  (4) 
IP(i,  0)  =  0,  Vi. 

The  operation  floor  (A)  denotes  the  nearest  integers  less  than  or 
equal  to  A.  The  results  of  the  projection  transform,  vectors  3ft  and 
©  will  have  the  dimension  of  nr  and  no,  respectively.  Both  projec¬ 
tions  3ft  and  ©  are  normalized  to  reduce  the  effect  of  illumination 
changes.  The  computation  of  projection  3ft  is  simple  and  does  not 
require  interpolation.  The  computation  of  projection  ©,  on  the  other 
hand,  requires  one-dimensional  interpolation,  as  shown  in  Eq.  (3). 

Examples  of  the  projections  of  APT  of  the  Lena  image  and  its 
scaled  and  rotated  version  are  shown  in  Fig.  3  (scale  =1.2  and  rota¬ 
tion  =45  degrees).  It  can  be  seen  that  the  scale  change  in  the  Carte¬ 
sian  appears  as  the  variable-scale  in  the  projection  3ft  (i.e.  from  f(t) 
to  f(at))  and  the  rotation  change  in  the  Cartesian  appears  as  shifting 
in  the  projection  ©. 

2.2.  PPT  Matching 

Given  a  model  image  IM(x,y),  a  target  image  IT(x,y)  that  is  a 
scaled  and  rotated  version  of  the  model  image  with  scale  parameter 
a  and  rotation  parameter  9  can  be  expressed  mathematically  as 

IT (x,y)  =  IM (ax  cos  9  +  ay  sin  9,—  ax  sin  9  +  ay  cos  9) .  (5) 
In  the  polar  domain,  Eq.  (5)  can  be  expressed  as: 

IT(r,e)  =  IM(ar,e  -6).  (6) 

Since  the  projections  3ft  and  ©  are  derived  from  PPT,  which  is  an 
approximation  of  Polar  transform,  we  can  obtain  scale  and  rotation 
parameters  between  two  images  by  comparing  their  projections.  As 
shown  in  Figs.  3(a)  and  3(c),  the  projections  ©  of  the  scaled  images 
in  the  Cartesian  are  slightly  altered  when  compared  with  that  of  the 
original  image.  This  is  because  the  areas  covered  during  the  PPT 
transformations  are  different  as  a  result  of  the  scaling.  Hence,  in 
order  to  accurately  obtain  the  shifting  parameter  between  the  two 
projections  ©M  and  ©T,  the  variable-scale  parameter  a  between  the 
two  projections  3 ftM  and  3ftT  needs  to  be  obtained  first. 
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Fig.  3.  The  effects  of  the  changes  in  scale  and  rotation  in  the  Carte¬ 
sian  coordinates  to  the  projections  9ft  and  ©:  (a)  the  scale  change  in 
the  Cartesian  appears  as  variable-scale  in  the  projection  9ft,  while  the 
projection  ©  becomes  slightly  altered,  (b)  the  rotation  in  the  Carte¬ 
sian  appears  as  shifting  in  the  projection  ©,  while  there  is  no  change 
in  the  projection  9ft,  and  (c)  the  changes  in  both  scale  and  rotation  in 
the  Cartesian  appear  as  variable-scale  in  projection  9ft  and  shifting  in 
projection  ©,  respectively 


2.2.1.  Find  scale  parameter 

There  are  several  ways  to  obtain  the  scale  parameter  from  the  pro¬ 
jections  depending  on  the  requirements  of  the  application  in  term 
of  the  computational  cost,  accuracy,  image  types  and  environments. 
We  introduce  here  two  effective  algorithms. 

•  Algorithm  1:  Logarithm  method 

This  algorithm  uses  the  scale  invariant  property  of  the  log¬ 
arithm  function.  First,  the  logarithm  function  is  applied  to 
the  projection  9ft  and  the  output  is  then  quantized  to  maintain 
the  original  dimension  of  the  projection.  The  mathematical 
expression  of  the  implementation  is  as  follow: 

OR(k)  =  5R(nr-i^SA);  k=l,...,nr.  (7) 

log  nr 

The  parameter  £9ft(.)  denotes  the  logarithmic  of  the  projec¬ 
tion  9ft.  Given  image  IT  sl  scaled  and  rotated  version  of  image 
/M,  the  scale  parameter  a  between  the  two  images  would  ap¬ 
pear  as  translation  in  the  logarithm  domain: 

mT(p)  =  mM  (j O )  +  log  a.  (8) 

To  find  the  displacement  d ,  where  d  =  log  a,  such  that 
£9ftT(p)  =  £9ftM(p  —  d),  one  can  evaluate  the  correla¬ 
tion  function  between  the  two  logarithmic  of  projections 

c(£jrm,  mTzy. 

d  =  argmaxC(£5RM,  £3tT).  (9) 

•  Algorithm  2:  Resampling  projection 

Algorithms  1  is  fast  but  can  yield  accuracy  at  some  degrees, 
since  the  displacement  estimation  can  be  computed  with  only 
integer  accuracy.  For  aerial  image  registration  purposes,  in 
which  high  precision  is  required,  one  may  include  resampling 
method  as  an  extra  step  after  acquiring  the  estimated  scale 


parameter  with  algorithm  1,  denoted  as  a! ,  (or  can  only  use 
resampling  method  if  estimated  scale  parameter  is  known  or 
expected  to  be  small).  Given  the  lower  limit  cll  and  the  up¬ 
per  limit  au  of  the  search  space  for  scale  parameter  such  that 
cll  <  a!  <  au,  respectively,  and  given  the  number  of  the 
search  steps  equals  to  h ,  we  resample  9ft^fs  =  T>(9 ftM,ai) 
where  a*  =  cll  +  t(^au~aL^  for  i  =  1,  ..h.  We  use  MSE 
method  to  compare  each  resampled  projections.  The  oper¬ 
ation  y  =  <F(x,  a)  is  a  resampling  procedure.  This  resam¬ 
pling  procedure  is  a  common  operation  in  signal  processing 
for  computing  the  sequence  in  vector  x  at  a  times  the  orig¬ 
inal  sampling  rate  by  using  a  polyphase  filter  implementa¬ 
tion.  Operation  T>  applies  an  anti-aliasing  linear-phase  (low- 
pass)  FIR  filter  to  x  during  the  resampling  process.  The  re¬ 
sult  will  have  the  dimension  of  vector  equal  to  length(y )  = 
length(x)  x  a.  The  scale  parameter  is  equal  to  the  resampling 
parameter  a%  that  yields  the  lowest  MSE,  Esr.  It  is  obvious 
that  the  larger  parameter  h  is,  the  higher  accuracy  would  be 
obtained. 

2.2.2.  Find  rotation  parameter 

After  the  scale  parameter  is  obtained,  the  next  step  is  to  find  the  rota¬ 
tion  parameter  0.  For  the  same  sampling  radius,  Rmax ,  the  larger  the 
image  (a  >>  1),  the  smaller  the  area  of  the  scene  is  covered  in  the 
sampling  procedure.  As  a  result,  the  magnitude  of  the  projections  0 
between  the  two  images  could  be  slightly  altered.  This  phenomenon 
is  illustrated  in  Figs.  3(a)  and  3(c).  Hence,  to  accurately  obtain  the 
rotation  parameter  0,  the  upper  limit  of  Eq.  (3)  needs  to  be  modi¬ 
fied  according  to  a.  If  a  >  1,  the  upper  limit  of  ©M  is  modified 
to  anr ,  while  the  upper  limit  of  computing  ©T  remains  unchanged. 
If  a  <  1,  the  upper  limit  of  ©T  is  modified  to  nr/a,  while  the  up¬ 
per  limit  of  computing  ©M  remains  unchanged.  Both  projections  are 
then  resampled  to  be  equal  in  length  before  comparison.  We  denoted 
the  modified  projections  as  ©M  and  ©T.  The  rotation  parameter  can 
be  found  by  evaluating  the  correlation  function: 

d  =  argmaxC(©M,  ©T),  6  =  2^.  (10) 

no 

2.2.3.  Find  distance  coefficient 

The  final  and  crucial  component  resulting  from  PPT  matching  is  the 
distance  coefficient,  denoted  as  S.  This  distance  coefficient  indicates 
how  large  the  difference  between  the  two  images  is.  The  distance 
coefficient  S  between  the  model  image  IM  and  the  target  image  IT 
can  be  computed  from  the  Euclidean  distance  between  the  projec- 

tions  0M  and  0T  as  £  =  ^E2i[©T(*)  -  -  d)]2. 

2.3.  Feature  Point  Based  Search  Scheme 

Similar  to  LPT  based  approach,  the  translation  parameter  between 
model  image  and  target  image  has  to  be  found  in  order  to  keep  the 
scale  and  rotation  invariant  properties.  We  adopt  Gabor  feature  ex¬ 
traction  to  accelerate  the  search  process  (more  detail  of  Gabor  fea¬ 
ture  extraction  can  be  found  in  author’s  publication  [4]).  These  fea¬ 
ture  points  are  obtained  by  applying  Gabor  transform  to  the  image 
and  selecting  those  pixels  which  have  high  energy  in  the  wavelet 
domain.  By  selecting  one  of  the  feature  points  in  the  model  image 
as  a  center  point  for  PPT,  we  reduce  the  search  from  every  pixel  of 
the  target  image  to  a  set  of  feature  points  found  in  the  target  image 
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only.  The  matched  feature  point  is  the  point  that  yields  the  lowest 
distance  coefficient  S.  Apparently,  the  number  of  feature  points  is 
much  smaller  than  that  of  the  pixels  in  the  target  image,  while  the 
computation  of  the  Gabor  wavelet  transform  is  much  lower  than  that 
of  PPT.  Thus,  the  computation  load  is  much  lighter  than  the  exhaus¬ 
tive  search  using  PPT. 

2.4.  Complete  algorithm 

Below  is  the  complete  steps  of  the  proposed  algorithm: 

•  Extract  feature  points  in  both  the  model  image  and  the  target 
image. 

•  Crop  circular  image  patch  IM  that  covers  the  area  Rmax  de¬ 
sired  to  be  registered  to  the  target  image. 

•  Compute  the  projections  9ftM  and  ©M  using  the  proposed 
PPT  approach  and  the  projection  transform. 

•  Use  each  feature  point  in  the  target  image  pj  for  z  =  1  :  nr, 
where  nr  is  the  number  of  feature  points  in  the  target  image 
as  the  origin  and  crop  a  circular  image  patches  if  for  z  =  1  : 
nr  with  the  radius  size  Rmax  • 

•  Compute  the  sets  of  candidate  projections  9ftT  =  {9ftf , 9ft  } 
and  ©T  =  {©f , ..,  @nT}- 

•  Match  each  candidate  with  the  model  using  the  proposed  PPT 
matching  algorithm.  Translation  is  the  point  (x,y)  that  yield 
the  lowest  distance  coefficient  S.  The  scale  and  rotation  pa¬ 
rameters  are  obtained  simultaneously  in  this  step. 

3.  EXPERIMENTAL  RESULTS 


(Scale  =1.01,  rotation  =13.89  degrees) 


(Scale  =1.00,  rotation  =-15.93  degrees) 


(Scale  =0.98,  rotation  =11.95  degrees) 


(a)  (b)  (c) 

(Scale  =1.02,  rotation  =49.95  degrees) 


Fig.  4.  Experimental  results  (a)  model  images,  (b)  target  images  and 
registration  results,  and  (c)  mosaic  images 


In  this  section,  the  performance  of  the  proposed  image  registration 
approach  is  evaluated.  We  uses  four  sets  of  aerial  images  in  the 
test  as  shown  in  Fig  4.  All  test  images  have  the  size  of  334  x  502 
pixels.  Since  the  scale  parameters  of  all  the  test  aerial  images  are 
very  small,  we  adopt  Algorithm  2  of  the  proposed  PPT  matching 
with  parameters  h  =  20,  cll  =  0.9  and  au  =  1.1  to  obtain  the  scale 
parameters.  The  size  of  the  search  space  using  Gabor  feature  point 
varies  depending  on  the  richness  of  details  in  the  image,  with  average 
of  120  -  150  points,  which  considers  as  approximately  0.1%  of  the 
total  number  of  pixels  in  the  image.  Figs.  4(a)  and  4(b)  represent  the 
model  image  and  the  target  image,  respectively.  In  Fig.  4(a)  the  red 
rectangular  area  indicates  the  image  patch  selected  to  be  registered 
to  the  target  image.  We  note  that  rectangular  shape  is  used  in  the 
figure  instead  of  circle  to  represent  the  orientation  of  the  images.  In 
Fig.  4(b),  the  red  rectangular  represents  the  registration  result.  To 
demonstrate  the  effectiveness  and  accuracy  of  the  proposed  method, 
we  create  mosaic  images,  in  which  the  model  images  are  scaled  and 
rotated  according  to  the  parameters  found  and  superimposed  onto 
the  target  images.  As  shown  in  Fig.  4(c),  one  can  see  that  images  are 
registered  accurately  and  naturally.  This  demonstrates  the  robustness 
to  translation,  scale  and  rotation  of  the  proposed  approach. 

4.  CONCLUSION 

In  this  paper,  we  have  presented  a  new  image  registration  approach 
using  the  projective  polar  transform.  With  the  combination  of  one- 
to-one  mapping  in  the  image  transform  process  and  the  innovative 
projection  approach  that  reduces  the  dimension  of  the  transformed 
image  from  two  to  one,  the  computational  load  in  both  mapping  and 
matching  process  is  very  low.  Two  effective  ID  matching  mech¬ 
anisms  are  proposed  to  effectively  and  accurately  obtain  scale  and 


rotation  parameters  of  the  registered  images.  The  experiments  with 

aerial  images  are  provided  to  demonstrate  the  accuracy  and  robust¬ 
ness  of  our  approach. 
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